Constraint damping in the Z4 formulation and harmonic gauge 
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We show that by adding suitable lower-order terms to the Z4 formulation of the Einstein equations, 
all constraint violations except constant modes are damped. This makes the Z4 formulation a 
particularly simple example of a A-system as suggested by Brodbeck et al. We also show that the 
Einstein equations in harmonic coordinates can be obtained from the Z4 formulation by a change 
of variables that leaves the implied constraint evolution system unchanged. Therefore the same 
method can be used to damp all constraints in the Einstein equations in harmonic gauge. 
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I. INTRODUCTION 

Formulations of the Einstein equations as an initial- 
boundary value problem are required for the numerical 
simulation of astrophysical events, such as the inspiral 
and merger of a binary system of black holes. Simula- 
tions are still crucially limited by instabilities. Most of 
these instabilities arise already in the continuum system 
of partial differential equations, rather that at the stage 
of finite differencing. 

There is now broad agreement that in order to avoid 
some of these instabilities, the formulation of the Einstein 
equations that one uses should give rise to well-posed 
problems. For the Cauchy problem (no boundaries in 
space) to be well-posed, strong hyperbolicity is a nec- 
essary and sufficient criterion. Symmetric hyperbolicity 
implies strong hyperbolicity and allows one to formulate 
a well-posed initial-boundary value problem. Symmet- 
ric hyperbolicity of the subsidiary constraint evolution 
system is likely to be a crucial ingredient in making the 
boundary conditions consistent with the constraints. 

In a hyperbolic formulation, the error associated with 
constraint violation grows at a bounded rate, but this 
can be very fast in practice. It would be preferable if 
one could find a formulation of the Einstein equations 
in which the submanifold of solutions that also obey the 
constraints is an attractor. Clearly this requires a mech- 
anism for breaking the time-reversal symmetry of gen- 
eral relativity away from the constraint surface. Mecha- 
nisms that have been suggested include dynamically ad- 
justing the free parameters of the constraint addition 0, 
or adding derivatives of the constraints so that the sys- 
tem becomes mixed parabolic and hyperbolic 0. 

Brodbeck et al 3] have suggested a general approach 
called the A-system to solving a system of evolution equa- 
tions and constraints such that the constraint surface is 
an attractor. This consists in adding one variable A for 
each constraint, such that the time derivatives of the A 
are the constraints, and the extended system is (or re- 
mains) symmetric hyperbolic. One then adds damping 
terms dtX = . . . — kX to the evolution equations for the 
variables A. As they are lower order, they do not af- 
fect the hyperbolicity of the system. These terms should 



damp the A, and therefore the constraints. 

In Q the Frittelli-Reula symmetric hyperbolic formu- 
lation of the Einstein equations was extended in this way, 
and it was shown analytically that when the system is 
linearised around Minkowski spacetime, the constraint 
surface is an attractor. In £| the conformal field equa- 
tions, reduced by two Killing vectors, were extended in 
the same way, and investigated numerically. This worked 
well for linear gravity, but not in strong field tests, where 
the constraint violations were reduced, but not to zero, 
and the error actually increased. 



II. Z4 AS A SIMPLE A-SYSTEM 

Our starting point is the observation that another sym- 
metric hyperbolic formulation of the Einstein equations, 
the Z4 formulation 0, 0, 13] is already a A system, with- 
out the need to add extra variables. (This was already 
noted in Q for the Z3 system, which is closely related 
to the Z4 system.) For simplicity we restrict our presen- 
tation to the vacuum case. Including matter would be 
straightforward. 

The Z4 system is obtained in its 4-dimensional co- 
variant form by replacing the vacuum Einstein equations 
Rab = by 



Rab + V a Z b + S7 b Z a = 0, 



(1) 



where R a b is the Ricci tensor of the 4-dimensional space- 
time metric g ab and Z a is an additional vector field. The 
main effect of this extension is to turn the 4 Einstein 
constraints into first-order evolution equations for the 4- 
vector Z a . A solution of the Z4 equations is a solution of 
the Einstein equations if and only if Z a = 0. (There is 
one exception: if the spacetime admits a Killing vector, a 
solution of the Z4 equations with Z a equal to the Killing 
vector is also a solution of the Einstein equations. In the 
following, we assume for simplicity that the spacetime is 
generic and does not admit a Killing vector.) 

We shall see that the variables Z a are already variables 
of the A type. All we need to do is to add the damping- 
term. We do this in covariant notation by replacing the 
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Einstein equations by 

R ab +V a Z b +V b Z a -n [t a Z b + t b Z a - (1 + P )g ab t c Z c ] = 0, 

(2) 

or, with the trace reversed, 

G ab +V a Z b +V b Z a -g ab V c Z c -K(t a Z b +t b Z a +pg ab t c Z c ) = 

where t a is a non- vanishing timelike vector field and n > 
and p are real constants. (We shall later restrict to 
p = 0.) It is £ a that explicitly breaks time reversal invari- 
ance. A simple geometrical choice is t a = n a , the future 
pointing unit normal vector on the time slicing. 
We carry out the usual 3+1 split of the metric as 

ds 2 = -a 2 dt 2 + 7 y (dx l + & dt) (dx j + ^ dt) , (4) 

and split Z a as Z a = X a + n a 9 where 9 = -n a Z a . In 
adapted coordinates this gives 8 = aZ° and Xi = Zi. 
(Note that n a as defined in is past-pointing, while 
ours is future-pointing, so that the two definitions of 9 are 
the same.) In the following we use the established, but 
slightly ambiguous, notation Z% to denote Xi. We denote 
the Ricci tensor of 7^ by Rij and its covariant derivative 
by £>,;. Spatial tensor indices i, j are moved with 7^. We 
also define the derivative operator do = a~ 1 (dt — Cp). In 
this notation, the 3+1 split of the Z4 equations is 

dolij = -2Kij, (5) 
d Kij = -a- l D i D j a + R ij -2K ik K k j +KK ij 

+DiZj + DjZi - 26Kij - «(1 + p)jije, (6) 

8 9 = X -H - 8K + D k Z k - D k (\na)Z k 

-{2 + p )k9, (7) 
d Zi = Mi + D t B - Di(lna)9 - 2K k Z k - nZ,. (8) 

In the last two equations H and Mi are shorthand for 
the Einstein constraints 

H = R - K i3 K lj + K 2 , (9) 
Mi = D j Kji - DiK. (10) 

What happens when the initial data for 7^ and Kij 
do not obey the Einstein constraints? By substituting 
the definitions of H and Mj into the evolution equations 
(0)-©, we find their formal time derivatives 

d Q H = -AM' 1 A (In a) + 2KB - 2D l M, - 4k(1 + p)K0 
+K ij ( 1 i ^ kl - 7 *V)[40l^fc - 4^fci], (11) 
doM, = --D i H+(K-20)Mi-Di(lna)H + 2R ij Z j 

+D 1 (D j Z l - DiZj) + Q-^ap^j + DjZi) 

-2D i (lna) J D i Z j - 2aT 1 K ii D 3 '[off) 

+2[K + K(l + p)]a- 1 D l (a9). (12) 

In order for the constraint evolution system to close (in 
the usual sense, see below) it must be considered to con- 
sist of <|11I12I7I8|) for the constraint variables H, Mi, 9 



and Zi. In particular, a solution of the evolution equa- 
tions obeys H = Mi = 9 = Z± = at all times if and 
only if they all vanish at t = 0. The constraint system 
associated with the Z4 system is unusual in that 9 and Zi 
are genuine dynamic variables while, as in other formula- 
tions, H and Mi are only shorthands for combinations of 
' the dynamic variables 7y and Kij and their derivatives. 
One can replace the 8 first-order evolution equations 
of the constraint system by a second-order wave equation 
for the 4- vector Z a by taking the divergence of 10, and 
using the contracted Bianchi identity W a G ab = 0. The 
result is 

OZ b + R ab Z a - nV a (t a Z b + t b Z a + pg ab t c Z c ) = 0, (13) 

where □ is the covariant wave operator V a V a . Given Z a , 
the Einstein constraints G 0ai can be read off from (J3J, or 
in 3+1 form H and Mi can be read off from 0-®- Note 
that for Zi = 8 = at some instant t = 0, the condition 
Zi — 8 = is equivalent to H = M, = at that instant. 
That means that all four constraints vanish at all times 
if and only if Zi = 9 = Zi = 9 = at t = 0. 

In either its first-order or second-order form, the con- 
straint evolution system closes only in the sense that the 
right-hand side is proportional to the constraints, but not 
in the sense that it is autonomous: one cannot consis- 
tently evolve either i|13fl or (|11I12I7I8|1 while considering 
the variables 7y and K^ (or equivalently the spacetime 
metric g ab ) as fixed. Instead one should focus on the 
following question: 

When one evolves initial data with a small constraint 
violation (set for example by finite differencing error) 
does the constraint violation grow or decay as the initial 
data are evolved? To address this, we perturb around a 
background solution gf b ' that obeys all 10 Einstein equa- 
tions, and write 

9a b = 9^+eg^, (14) 

where IT? = and Z^ =0. To first order in e the 
constraint violation then obeys a linear evolution equa- 
tion with coefficients taken from the background Einstein 
spacetime, and admitting arbitrary data. In the Z4 for- 
mulation this is just a vector wave equation on the back- 
ground spacetime, namely 

□ (°)Z«- K V«W (t^ZP +pg^t^Z^) = 0. 

(15) 

The growth of sufficiently small constraint violations is 
controlled by the damping; for larger constraint viola- 
tions nonlinear lower-order terms also become important. 
We expect that with sufficiently large n, and starting 
from sufficiently small initial constraint violations, the 
nonlinear terms never become important in practice. 

[Note added in revision: Friedrich has has indepen- 
dently carried out a partial non-linear analysis of l|13|) 
in the context of generalised harmonic coordinates and 
without damping. In 113|l he considers the spacetime in 
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the wave operator □ as given, but expresses R a b in terms 
of Z a using @ ■ The result is wave equation for Z a with a 
quadratic lower-order term, but on a fixed spacetime. He 
shows that generic solutions blow up in finite time. This 
analysis accounts for some of the hidden nonlinearity of 
(|13|1 . but does not include all terms of 0(e 2 ). We have 
restricted ourselves to a consistent linear analysis.] 

III. MODE ANALYSIS 

The constraint evolution system of the Z4 system is 
simpler than in other formulations of the Einstein equa- 
tions in that it has the form of a covariant wave equation. 
We can use this to carry out a mode analysis for the lin- 
earised constraint system (I15fl on an arbitrary Einstein 
background in the high-frequency, frozen coefficient limit 
in which the wavelength of Z a is much smaller than the 



curvature scale of the Einstein background. We can then 
locally approximate the background g$ as Minkowski 
space, and work in standard Minkowski coordinates. As- 
suming without loss of generality that t a t a — —1, we 
go to the frame in which i M = (1,0,0,0) and find (now 
dropping the expansion indices) 

OZ - k[(2 + p)d t Z -d i Z i ] = 0, (16) 
nZi-K(dtZi + pdiZ ) = 0, (17) 

where □ is now the Minkowski wave operator —df + did 1 . 
We make a plane-wave ansatz 

Z^t,x i ) = e st+ ^ i Z^, (18) 

with complex s and real u>i. This gives rise to the eigen- 
value problem 



2 - k(2 + p)s niuj 

-npiuu -s 2 -uj 2 - ks ) | Z n ) = 0. (19) 

-S 2 — Ld 2 — K 




where Z n is the component of Zi in the direction of u>i 
and Zj± stands for the projection of Zi normal to w,. The 
4 eigenvalues s for Za are 

■ = -5*1/(5)'-"" (20) 

(each of these occurring twice), independently of p. The 
other 4 eigenvalues s are in general complicated. How- 
ever, for p = they take the simple form 

s = -h±\/k 2 -uj 2 (21) 

(each of these occurring twice) . Therefore with p = all 
modes are damped for all uji ^ 0. At high wave numbers, 
u) » k, the damping is by a constant factor, with 

K 

s ~ — k ± iuj, — — ± iu>, (22) 

but at low wave numbers, uj -C k, it is similar to a heat 
equation, with 

uj 2 uj 2 
«~-«,— -,-2/s,-— . (23) 

Half of the modes are damped less with decreasing 
wavenumber, and not at all at zero wavenumber. The 
only other case in which the eigenvalues s are simple is 
p = — 1. This also simplifies the field equations, but in 
this case the other 4 eigenvalues are s — —n±iuj and 
s = ±iuj, and so there are undamped modes for all Wj. 



It turns out that all but the constant modes are damped 
for any p > — 1, but p = is the most natural choice, 
and we assume this value in the following. 

We have shown that the constraint manifold is an 
attractor for k > when the equations are linearised 
around Minkowski space, with the exception of constraint 
violations that are constant in space. The same is true 
in the high-frequency limit when linearising around any 
Einstein background. This analysis breaks down when 
the constraint violations are large and/or when their 
wavelength is large compared to the background curva- 
ture. In that case lower-order terms can potentially make 
the constraints grow against the explicit damping term. 
However, compared to the systems of yJ-LI- ZA has only 
a fraction of the number of variables and constraints, and 
so is less likely to be affected by undesirable growth aris- 
ing from lower order terms. 



IV. Z4 AND HARMONIC GAUGE 

The hyperbolicity of any formulation of the Einstein 
equation depends on the choice of gauge (assuming that 
the formulation does not already fix the gauge). Formu- 
lations of the Einstein equations derived from the ADM 
formulation are typically not hyperbolic when the lapse 
a and shift j3 l are given functions of the coordinates, but 
if they can be made hyperbolic at all, this is typically 
true for fixed shift and fixed densitised lapse Q, where a 
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and Q are related by 

a = j^ 2 Q, (24) 

for some constant parameter a. The Z4 formulation with 
fixed densitised lapse and fixed shift is strongly hyper- 
bolic for a > 0. Surprisingly, it is not symmetric hyper- 
bolic for any a. (The most general energy that is con- 
served in the high-frequency, frozen coefficient approxi- 
mation fails to be positive definite. Details will be given 
elsewhere.) 

An evolved version of the densitised lapse (usually 
called the Bona-Masso lapse) is 

d In a = -a(K - mO) + d In Q, (25) 

where m is another constant parameter. When 9 = 0, 
the 8q derivative of l|24|l is just i|25|) . As pointed out in 
0, the Z4 formulation is symmetric hyperbolic with the 
lapse (I25[) and fixed shift for a = 1 and m = 2. For any 
other a > 0, it is strongly hyperbolic (with arbitrary m), 
but not symmetric hyperbolic. (The most general energy 
fails to be positive definite for a ^ 1. Details will be given 
elsewhere.) a = 1 is equivalent to harmonic slicing. This 
suggests that in some way Z4 is closely related to the 
Einstein equations in harmonic coordinates. 

In fact, if in we consider Z^ as the shorthand 

Z» = \ (-H, - 5 a/3 W + \g ap g a p^ ■ (26) 

where are given functions of the coordinates ( "gauge 
source functions"), we obtain the Einstein equations in 
the generalised harmonic gauge. In harmonic gauge all 
ten g^ v , or equivalently 7^, a and (3 l obey quasilinear 
second-order evolution equations whose principal part is 
just the wave operator on the spacetime with metric g^ u . 

The constraints = must still be obeyed to obtain 
a solution of the Einstein equations, but they are now 
the harmonic gauge constraints 

- <r >v + \g^g aP g aP ,„ = H», (27) 

where □ is the scalar wave operator [T(|. In a 3+1 split, 
these constraints can be solved for a and $ l , and thus 
constrain the free data for the wave equations for a and 

The substitution (|26|l does not change the constraint 
system itself at all, which is still the wave equation for 
Z^. The only difference is now that the Z^ are no longer 
dynamical variables but are shorthands for the harmonic 
gauge constraints. We immediately obtain a prescription 
for damping all constraints (gauge and Einstein) in nu- 
merical free harmonic evolutions: we modify the Einstein 
equations in harmonic gauge to 

Rpv + V M Z„ + V U Z M - n{t^Z v + tyZ^ - g^ut Z\) = 0, 

(28) 

where Z^ is the shorthand (|26|l . (We have used Greek in- 
dices to indicate that these are not tensor equations but 



hold only in harmonic coordinates.) As the constraint 
system is identical to that of Z4, the same mode analy- 
sis applies, showing that all Fourier modes are damped 
except the constant in space modes. 

V. Z4 AND NOR/BSSN 

We obtain the NOR system by the substitutions 

0-0, Z,-^(/,-1 J V + ^7 Al ) (29) 

in terms of a new variable fj. (p is the constant parameter 
of the same name in |llj . but is different from the p 
introduced in (J2J, which is now set to zero.) The variable 
9 disappears, and Zi is now the definition constraint of 
the new auxiliar y v ariable L . (Note that Zi = —Gi in 
the notation of [ll| andllj-) Essentially this change 
of variables was used in [fj to obtain the BSSN system 
from the Z4 system, where constraint damping was also 
pointed out. The constraints are H = 0, Mj = 0, and 
Zi = 0. The constraint evolution system is obtained 
from that of the Z4 system by setting 9 to zero. (The Z4 
system with 9 = but treating Zi as dynamical variables 
is called the Z3 system, and was proposed in 0. Its 
constraint evolution system is identical to that of the 
NOR system.) 

If we repeat the mode analysis for the Z3 system, we 
find that the dispersion relation relation for the group of 
vectors transverse to u>i, (Ma, Z a) is given by 

s = ~(-k± V* 2 - 4w 2 ) , (30) 

but for the group of scalars (H, M n , Z„) the modes are 
given by 

s = — k, ±iu>, (31) 

so there are two undamped modes. This means we can 
damp 5 of the 7 constraints of the NOR system (or the 
BSSN system) by a suitable modification given essentially 
by (H with 113 • 

VI. CONCLUSIONS 

We have given a A system for general relativity that is 
optimal in the sense that the only variables in addition 
to the usual ADM variables 7^ and Kij are the four A 
variables 9 and Zi that are associated with the four Ein- 
stein constraints H and Mj. All constraints except the 
constant in space modes are damped through lower or- 
der terms. This system is already symmetric hyperbolic 
(for harmonic slicing and fixed shift) without the need for 
other auxiliary variables. It is therefore a good testbed 
for investigating the usefulness of A systems in general 
relativity in general. 
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Within the framework of any A system, it is impossi- 
ble to obtain damping of homogeneous constraint modes 
because the damping is through a friction term. How- 
ever, it is likely that in the near future well-posed initial- 
boundary value problems will be constructed in which 
the constraint system implicitly obeys maximally dissi- 
pative boundary conditions. When these are of Dirichlet 
type, the homogeneous constraint modes should also be 
eliminated. 

We have pointed out that the Z4 system is closely re- 
lated to harmonic gauge, and that the two share the same 
constraint evolution system. Therefore constraint damp- 
ing works in the same way in harmonic gauge evolutions, 



and all 8 constraints can be damped. Similarly, in the Z3 
system and its counterpart the NOR/BSSN system, 5 of 
7 constraints can be damped. 
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